运筹系列51：python |
您所在的位置：网站首页 › python 整数规划 › 运筹系列51：python |
运筹系列51：python
|
1. 介绍
CBC和Gurobi的一个python封装，支持cut generation, lazy constraints, MIP starts以及solution pools。 自带CBC，省去安装的麻烦。 python-mip使用cffi库调用C代码，使用Pypy进行编译，据称性能很好，比gurobi自带的python接口快25倍，比JuMP还要快（但总体上差的不多）。
定义模型 m = Model() # default is CBC m = Model(sense=MAXIMIZE, solver_name=GRB) # use GRB for Gurobi添加目标函数： m.objective = xsum(c[i]*x[i] for i in range(n)) # 默认是min m.objective = maximize(xsum(c[i]*x[i] for i in range(n))) # 手动添加添加变量： x = m.add_var() # 默认是continues，此外还有binary和integer。此外还可以定义lb和ub，以及给变量命名 z = m.add_var(name='zCost', var_type=INTEGER, lb=-10, ub=10) # 默认是非负的，除非人工设置lb y = [ m.add_var(var_type=BINARY) for i in range(10) ] # 变量数组 vz = m.var_by_name('zCost') # 获取变量的方法 vz.ub = 5添加约束，非常独特，没有方法，直接用+ m += x + y = y[j]-n # optimizing model.optimize() # checking if a solution was found if model.num_solutions: out.write('route with total distance %g found: %s' % (model.objective_value, places[0])) nc = 0 while True: nc = [i for i in V if x[nc][i].x >= 0.99][0] out.write(' -> %s' % places[nc]) if nc == 0: break out.write('\n') 3.3 location：引入SOS变量 import matplotlib.pyplot as plt from math import sqrt, log from itertools import product from mip import Model, xsum, minimize, OptimizationStatus # possible plants F = [1, 2, 3, 4, 5, 6] # possible plant installation positions pf = {1: (1, 38), 2: (31, 40), 3: (23, 59), 4: (76, 51), 5: (93, 51), 6: (63, 74)} # maximum plant capacity c = {1: 1955, 2: 1932, 3: 1987, 4: 1823, 5: 1718, 6: 1742} # clients C = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] # position of clients pc = {1: (94, 10), 2: (57, 26), 3: (74, 44), 4: (27, 51), 5: (78, 30), 6: (23, 30), 7: (20, 72), 8: (3, 27), 9: (5, 39), 10: (51, 1)} # demands d = {1: 302, 2: 273, 3: 275, 4: 266, 5: 287, 6: 296, 7: 297, 8: 310, 9: 302, 10: 309} # plotting possible plant locations for i, p in pf.items(): plt.scatter((p[0]), (p[1]), marker="^", color="purple", s=50) plt.text((p[0]), (p[1]), "$f_%d$" % i) # plotting location of clients for i, p in pc.items(): plt.scatter((p[0]), (p[1]), marker="o", color="black", s=15) plt.text((p[0]), (p[1]), "$c_{%d}$" % i) plt.text((20), (78), "Region 1") plt.text((70), (78), "Region 2") plt.plot((50, 50), (0, 80)) dist = {(f, c): round(sqrt((pf[f][0] - pc[c][0]) ** 2 + (pf[f][1] - pc[c][1]) ** 2), 1) for (f, c) in product(F, C) } m = Model() z = {i: m.add_var(ub=c[i]) for i in F} # plant capacity # Type 1 SOS: only one plant per region for r in [0, 1]: # set of plants in region r Fr = [i for i in F if r * 50 = 1e-6]: plt.plot( (pf[i][0], pc[j][0]), (pf[i][1], pc[j][1]), linestyle="--", color="darkgray" ) plt.savefig("location-sol.pdf") 4. 自定义约束动机：放松约束容易求解。每次迭代时需要找到最违背的约束，这也是个优化问题，叫做separation problem。 4.1 Cutting plane，也叫行生成TSP问题，subtour elimination用的如下形式： 接下来的问题是，每次如何获得这些约束？使用 Minimum cut problem in graphs，在python中可以使用networkx min-cut module from itertools import product from networkx import minimum_cut, DiGraph from mip import Model, xsum, BINARY, OptimizationStatus, CutType N = ["a", "b", "c", "d", "e", "f", "g"] A = { ("a", "d"): 56, ("d", "a"): 67, ("a", "b"): 49, ("b", "a"): 50, ("f", "c"): 35, ("g", "b"): 35, ("g", "b"): 35, ("b", "g"): 25, ("a", "c"): 80, ("c", "a"): 99, ("e", "f"): 20, ("f", "e"): 20, ("g", "e"): 38, ("e", "g"): 49, ("g", "f"): 37, ("f", "g"): 32, ("b", "e"): 21, ("e", "b"): 30, ("a", "g"): 47, ("g", "a"): 68, ("d", "c"): 37, ("c", "d"): 52, ("d", "e"): 15, ("e", "d"): 20, ("d", "b"): 39, ("b", "d"): 37, ("c", "f"): 35, } Aout = {n: [a for a in A if a[0] == n] for n in N} Ain = {n: [a for a in A if a[1] == n] for n in N} m = Model() x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A} m.objective = xsum(c * x[a] for a, c in A.items()) for n in N: m += xsum(x[a] for a in Aout[n]) == 1, "out({})".format(n) m += xsum(x[a] for a in Ain[n]) == 1, "in({})".format(n) newConstraints = True while newConstraints: m.optimize(relax=True) print("status: {} objective value : {}".format(m.status, m.objective_value)) G = DiGraph() for a in A: G.add_edge(a[0], a[1], capacity=x[a].x) newConstraints = False for (n1, n2) in [(i, j) for (i, j) in product(N, N) if i != j]: cut_value, (S, NS) = minimum_cut(G, n1, n2) if cut_value |
距离矩阵如下： 
将第三个约束分离（separate），得到如下解： 
添加两条约束：𝑥(𝑑,𝑒)+𝑥(𝑒,𝑑)≤1 以及 𝑥(𝑐,𝑓)+𝑥(𝑓,𝑐)≤1，求解得到如下 
增加约束：𝑥(𝑎,𝑔)+𝑥(𝑔,𝑎)+𝑥(𝑎,𝑏)+𝑥(𝑏,𝑎)+𝑥(𝑏,𝑔)+𝑥(𝑔,𝑏)≤2 求解得到结果为261。【本文地址】
今日新闻 |
推荐新闻 |